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ABSTRACT 

We investigate how well the intrinsic shape of early-type galaxies can be recovered when 
both photometric and two-dimensional stellar kinematic observations are available. We simulate 
these observations with galaxy models that are representative of observed oblate fast-rotator 
to triaxial slow-rotator early-type galaxies. By fitting realistic triaxial dynamical models to 
these simulated observations, we recover the intrinsic shape (and mass-to-light ratio), without 
making additional (ad-hoc) assumptions on the orientation. 

For (near) axisymmetric galaxies the dynamical modelling can strongly exclude triaxiality, 
but the regular kinematics do not further tighten the constraint on the intrinsic flattening 
significantly, so that the inclination is nearly unconstrained above the photometric lower limit 
even with two-dimensional stellar kinematics. Triaxial galaxies can have additional complexity 
in both the observed photometry and kinematics, such as twists and (central) kinematically 
decoupled components, which allows the intrinsic shape to be accurately recovered. For 
galaxies that are very round or show no significant rotation, recovery of the shape is degenerate, 
unless additional constraints such as from a thin disk are available. 

Key words: galaxies: elliptical and lenticular, cD - galaxies: kinematics and dynamics - 
galaxies: structure 



1 INTRODUCTION 

Numerical cold dark matter simulations predict that the dark mat- 
ter haloes of galaxies typically have intrinsic shapes with three 
distinct axes a, b and c (e.g. Jing & Suto 2002). Such a triaxial 
shape is characterized by the axis ratios p = b/a and q = c/a and 
^ q p ^ 1, which can vary as a function of radius. Based on 
the specific angular momentum inferred from integral-field stellar 
kinematics, early-type galaxies seem to be divided into two classes 
of fast and slow rotators, different from their morphological clas- 
sification in ellipticals and lenticulars (Kormendy & Bender 1996; 
Emsellem et al. 2007). Most fast rotators, including lenticular as 
well as many elliptical galaxies, which show strong rotation around 
their photometric minor axis, seem consistent with oblate axisym- 
metry (p = 1). However, many of these systems contain bars and 
also weakly triaxial shapes can not be ruled out. On the other hand, 
the slow rotators, with little or no rotation, are most likely triaxial. 

However, all this is based on observed quantities, i.e., projected 
onto the plane of the sky, while for the true comparison between 
galaxies intrinsic properties, and thus the viewing directions are 
needed. Most early-type galaxies lack thin discs (in either dust, gas 
or stars), so that typically there is only a weak constraint on the 
viewing direction from (mainly) the observed ellipticity. 

* Hubble Fellow 



Determining the intrinsic shape of galaxies is possible by statis- 
tical analysis of their observed ellipticities (e.g. Lambas, Maddox & 
Loveday 1992; Ryden 1992). The deficiency of E0 type galaxies im- 
mediately rules out that elliptical are perfectly axisymmetric (Franx, 
Illingworth & de Zeeuw 1991; Ryden 1996), even if bars are taken 
into account. A recent study by Padilla & Strauss (2008) using ellip- 
ticities from SDSS show that early-type galaxies are weakly triaxial 
with mean axial ratios of p ~ 0.95 and q ~ 0.7. However, these 
kind of analyses involve rather strong assumptions on the intrinsic 
shape distribution, and can not determine the shape of individual 
galaxies. Alternatively, by assuming a functional form for the intrin- 
sic density it is possible to constrain the intrinsic shape of the light 
distribution of an individual galaxy, given its observed ellipticity 
and isophotal twist profile (e.g. Williams 1981; Chakraborty et al. 
2008). In general, however, there is a well-known non-uniqueness 
in the deprojection of the surface brightness — even for axisym- 
metric galaxies not viewed perpendicular to their symmetry axis — 
also called the 'cone of ignorance' (Rybicki 1987). This tells us 
that it is not possible to constrain the intrinsic shape using only the 
photometry without making strong assumptions. 

Nowadays, in addition to observations of the surface bright- 
ness of early-type galaxies, integral-field spectrographs provide the 
full line-of-sight velocity distribution (LOSVD) as a function of 
position on the sky C(x' , y' , v z i). The LOSVD is commonly pa- 
rameterized in terms of its velocity moments, resulting in maps of 
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the mean line-of-sight velocity, velocity dispersion and higher-order 
velocity moments. Can we (better) constrain the intrinsic shape, and 
hence the viewing direction, of galaxies from these two-dimensional 
photometric and kinematic observations? 

Because the deprojection is degenerate and the conversion from 
light to mass is complicated when dark matter is present, imply- 
ing a range of possible gravitational potentials. Even ignoring (or 
parameterizing) both these aspects, it is far from evident that the 
LOSVD provides enough information to infer the intrinsic shape. 
Let us first consider three-integral oblate axisymmetric models, i.e., 
with a distribution function (DF) that is a function f(E, L z ,h) of 
the energy E, the angular momentum component L z parallel to the 
symmetry z-axis, and a non-classical third integral of motion 1$. In 
the special case that a galaxy happens to be well approximated by 
a two-integral DF f(E, L z ), the density p(R, z) (in the cylindrical 
coordinates R and z) uniquely determines the even part of f(E, L z ) 
and the mean streaming p{v$) in the meridional plane fixes the part 
of f(E, L z ) that is odd in L z (Dejonghe 1986). If both p(R, z) and 
p{vj,) are known, they define a two-integral DF completely. The 
observed velocity dispersion and higher order velocity moments 
of the LOSVD then provide additional information, which for ex- 
ample can be used to constrain the inclination (see e.g. Cappellari 
et al. 2006). This two-integral semi-isotropic result was recently 
generalized by Cappellari (2008) using a special class of anisotropic 
models in which the velocity ellipsoid is assumed to be aligned with 
cylindrical coordinates, and the anisotropy f3 z (as defined by Binney 
& Mamon 1982) to be positive and constant within a galaxy. Both 
assumption are based on an simplified interpretation of fig. 1 in 
Cappellari et al. (2006), and is at most a first-order approximation. 

In the realistic case of a three-integral DF f(E,L z ,Iz), it 
might well be that the full three-variable LOSVD C{x' , y' ,v z i) is 
needed to determine the DF, so that there is no information left to 
constrain the inclination. Using long-slit spectroscopy it was gen- 
erally impossible to constrain the inclination (e.g. van der Marel 
et al. 1998; Gebhardt et al. 2000). Initially, the availability of two- 
dimensional stellar kinematics seem to allow a constraint on the 
inclination beyond the lower limit from the photometry, even in the 
case of full three-integral axisymmetric orbit-based Schwarzschild 
models (Verolme et al. 2002). However, Krajnovic et al. (2005) con- 
vincingly showed that, after carefully taking into account numerical 
effects, similar three-integral axisymmetric Schwarzschild models 
at different inclinations above the photometric limit can (nearly) all 
fit the observed LOSVD within the measurement uncertainties. 

Thomas et al. (2007) showed that axisymmetric modeling can 
have signficant systematics in the recovered mass-to-light ratios 
by modelling triaxial merger remnants with their axisymmetric 
machinery. This shows that triaxial modelling is strongly preferred 
over axisymmetric modelling, as real elliptical galaxies are at least 
weakly triaxial. 

In the triaxial case, the DF is again a function of three integrals 
of motion f(E, I2, h), but the orbital structure in these models is 
substantially richer than in the oblate axisymmetric models, with 
four major orbit families, instead of only one. This introduces a 
"fundamental" non-uniqueness in the recovery of the DF: whereas in 
the oblate axisymmetric case p(R, z) uniquely defines the even part 
of f(E, L z ), in the (separable) triaxial case the density p(x, y, z) 
does not uniquely determine the even part of f(E, h, even 
though both of these are functions of three variables (Hunter & 
de Zeeuw 1992). It is (yet) unknown how much specification of 
C(x' , y ,v z i) can narrow down the range of possible DFs further. 
Therefore, in the triaxial case it seems to be even harder, if not 
impossible, to constrain the viewing direction. 



However, while mathematically there are an infinite number of 
deprojections of the surface brightness possible already for oblate 
axisymmetric systems viewed away from edge-on, in practice only a 
limited range of deprojections lead to intrinsic density distributions 
that are considered realistic for galaxies. Similarly, we have shown 
in (van de Ven, de Zeeuw & van den Bosch 2008, hereafter vdV08) 
that the DF in the triaxial case is in practice well recovered at the 
correct viewing direction using orbit based models. At the same 
time, we found in (van den Bosch et al. 2008, hereafter vdB08) that 
if there is enough complexity in the photometry and/or kinematics 
it is likely that even the intrinsic triaxial shape is well constrained. 
Such complexities are for example photometric and kinematic twists, 
misalignment between the photometric minor axis and kinematic 
rotation axis (kinematic misalignment) and so-called kinematically 
decoupled components (KDCs). For example, the velocity field of 
NGC4365 (Davies et al. 2001) shows rotation around the photo- 
metric short-axis in the core as well as around the long-axis in the 
outer parts. While the central KDC excludes a prolate shape, the 
global long-axis rotation is not possible in an oblate system, leaving 
a triaxial intrinsic shape as the only solution (see also vdB08). 

Statistical analysis of observed ellipticities and kinematic mis- 
alignments can provide an estimate of the intrinsic shape distribution 
(Binney 1985; Franx et al. 1991). However, this statistical approach 
requires a large sample of galaxies for which the kinematic misalign- 
ment is accurately measured, which is not (yet) available. Tenjes et al. 
(1993) explored the use of the observed ellipticity and kinematic 
misalignment to constrain the intrinsic shape of individual galax- 
ies, adopting a specfic form for the intrinsic density and streaming 
motion. Statler (1994, and references therein) developed a velocity- 
field fitting method for individual galaxies based on a solution of 
the continuity equation, assuming a similar flow which separates 
the density and streaming motion in a radial and angular part. Be- 
cause this approach is fast the large parameter space can be sampled 
finely enough to find the most probable solution using Bayesian 
statistics. However, the method assumes a (plausible) solution for 
the continuity equation, is not self-consistent, and fits the first ve- 
locity moment only. By comparison, the orbit-based Schwarzschild 
models that we employ allow for general density distributions, are 
self-consistent, and fit the full LOSVD. As mentioned above, even 
with such a general method making optimal use of the state-of-the- 
art two-dimensional photometric and kinematic observations, it is 
not evident how well the intrinsic shape can be recovered. 

In this paper, we use a set of realistic (nearly) axisymmetric 
and triaxial galaxy models with a three-integral DF, to show that 
an increasing complexity in the corresponding two-dimensional 
observables indeed leads to a better recovery of the intrinsic shape. 
In § 2, we construct these Abel models based on the description 
of vdV08. We fit the corresponding observables using the triaxial 
Schwarzschild method of vdB08, of which we highlight the most 
relevant aspects in § 3. We investigate the resulting constraints on 
the intrinsic shape in § 4, and briefly discuss the recovery of the 
intrinsic moments in § 5. We discuss our findings and conclude in 
§6. 



2 GALAXY MODELS WITH A THREE-INTEGRAL 
DISTRIBUTION FUNCTION 

Following vdV08, we build triaxial Abel models with a separable (or 
Stackel) potential that are a generalization of the spherical Osipkov- 
Merritt (Osipkov 1979; Merritt 1985) models (see also Dejonghe 
& Laurent 1991; Mathieu & Dejonghe 1999). These models, with 
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a DF f(E, 12,13) that depends on three exact integrals of motion, 
allow for a large range of shapes and internal dynamics, while the 
corresponding observables, including the LOSVD, can be computed 
efficiently. They are well suited to simulate realistic imaging and 
integral-field kinematics of early-type galaxies. Here, we construct 
thirteen galaxy models, six nearly oblate axisymmetric, three nearly 
prolate axisymmetric and four triaxial. These galaxy models are 
representative of early-type galaxies observed with the integral-field 
spectrograph SAURON (e.g. Emsellem et al. 2004), from oblate 
fast-rotator to triaxial slow-rotator. 

As in § 4 of vdV08, we use a triaxial isochrone Stackel potential 
with a length scale of 10" ~ 1 kpc at the adopted distance of 
20Mpc, and with a total mass of 10 11 Mq. Next, we choose the 
DF to be a linear combination of components of the form f(S) = 
(S - 5iim)/(l - Sim,), where both S = -E + wl 2 + ul 3 and 
Sum depend on the (constant) shape parameters w and u. Moreover, 
each DF component can be one of three types. The non-rotating 
(NR) type is made of box orbits and tube orbits with both senses of 
rotation populated equally. In addition we have two rotating types, 
long-axis rotating (LR) and short-axis rotating (SR), which consist 
of tube orbits, and have net rotation around the long axis and short 
axis respectively. We set u and w of each component such that the 
total density is similar to that corresponding to the potential, except 
in the outer parts (> 30") where a steeper surface brightness profile 
mimics the presence of dark matter (see also vdV08). 

To calculate the observables we need to choose a viewing 
direction. For the oblate/prolate axisymmetric models we adopt 
three inclinations: close to edge-on/side-on i = 87° (cosi ~ 0.05), 
intermediate i = 60° (cosi = 0.50), and close to face-on/end-on 
i — 18° (cosi ~ 0.95). For the triaxial models we use polar and 
azimuthal viewing angles of $ = tp = 60°. We then compute 
for each DF component on a rectangular grid on the sky plane the 
LOSVD and fit a Gauss-Hermite series to obtain maps of the surface 
mass density S, and the mean line-of-sight velocity V, velocity 
dispersion a and higher-order Gauss-Hermite velocity moments up 
to he (Gerhard 1993; van der Marel & Franx 1993). The surface 
brightness is obtained from dividing E by a constant (total) mass-to- 
light ratio M/L= 4 M /L . To simulate SAURON observations, we 
combine pixels to obtain a minimum signal-to-noise (S/N), which 
we take proportional to the square root of the surface brightness. 
Adopting a typical mean error of 7.5 km s _1 in V and <r, and 0.03 
for hs through h,6, we assign per (Voronoi) bin an error, weighted 
with the S/N, and use these errors to (Gaussian) randomize the 
kinematic maps. 

The resulting maps of the surface brightness, V, a, hz and hi 
are shown in top rows in Figures 1-5 for each of the thirteen galaxy 
models, for which the specific choices are as follows. 

Fast rotator oblate models (F01,F02,F03) We choose the 
isochrone Stackel potential so that the (mean) axis ratios of the 
corresponding density are ps = 0.99 and qs = 0.76, such that 
the model is close to oblate axisymmetric 1 . The DF contains a 
NR and SR component with relative mass fractions of 30 and 70 
per cent, and both with shape parameters w = u = —0.5. This 
results in regular velocity field around the (photometric) minor 
axis. From a luminosity weighted average of the simulated V and 
a maps, we calculate the (projected) specific angular momentum 
Xr = (R\V\)/(R\/V* +cr 2 ) (see Emsellem et al. 2007). Both the 
edge-on and intermediate inclination models FOl and F02, have 



1 In the axisymmetric limit the isochrone potential reduces to the Kuzmin 
& Kutuzov (1962) potential (cf. Dejonghe & de Zeeuw 1988). 



Xr values of 0.36 and 0.33 that are well within the regime of the 
fast rotators. Even the close to face-on model F03 has a Xr value of 
0.14 that is above the boundary Xr = 0.10 between fast and slow 
rotators. 

Slow rotator oblate models ( S01,S02,S03) These three models 
are equivalent to the above fast rotator oblate models, except that 
the relative mass fractions of the DF components are inverted, i.e., 
70 and 30 per cent for the NR and SR component, respectively. As a 
result, they rotate significant slower, with Xr values for the edge-on, 
intermediate and face-on inclination models SOI, S02 and S03 of 
respectively 0.16, 0.15 and 0.06. 

Slow rotator prolate models (SP1,SP2,SP3) The isochrone 
potential in this case is chosen such that the axis ratios of the cor- 
responding density are ps = 0.79 and qs = 0.78, i.e., close to 
prolate axisymmetric. The DF consists of a NR and LR component 
with relative mass fractions of 70 and 30 per cent, again both with 
shape parameters w = u = —0.5. The resulting kinematics show 
regular rotation around the (photometric) major axis. The amplitude 
of rotation is rather mild, with A r values for the side-on, intermedi- 
ate and end-on models SP1, SP2 and SP3 of respectively 0.18, 0.15 
and 0.05. 

Slow rotator triaxial models (ST1 ,ST2,ST3,ST4) The axis ratios 
ps = 0.90 and qs = 0.77 of the density corresponding to the 
chosen isochrone potential are consistent with a triaxial model. 
Model ST1 only has a NR component with w = u = —0.5, and 
hence shows no rotation (Xr = 0), similar to NGC 4486 (M 87). 
Model ST2 has an additional SR component with relative mass 
fraction of 20 per cent, so that X R — 0.14. Model ST3 has in 
addition to the NR component a SR component with relative mass 
fraction of 10 per cent and different shape parameters w — 0.5 and 
u = —1.0. This results in a central compact kpc-size KDC and zero 
rotation outwards due to the extended NR component, similar to 
slow rotators like NGC 5831. The Xr value within the inner 10" 
rises above 0.10, but at larger radii quickly drops within the slow 
rotator regime. Finally, model ST4 is the same as model ST3, but 
with an additional LR component with u — w — —0.5. The three 
NR, SR and LR components have relative mass fractions of 80, 
10 and 10 per cent. Outside the KDC, the model has significant 
(long-axis) rotation, similar to NGC 4365. As a result, the Xr value 
stays around 0.10 in the outer parts. 

Because the Abel models are constructed with a minimum num- 
ber of DF components, they are transparent and relatively quick to 
calculate, but still they capture most of the rich dynamics observed 
in early-type galaxies. Moreover, a careful decomposition of the 
observed SAURON two-dimensional stellar kinematics of early-type 
galaxies using kinemetry shows that most early-type galaxies can be 
described by only a few components (Krajnovic et al. 2008). Still, 
due to the sometimes sharp transition between the DF components, 
the surface brightness of the Abel models is not always as smooth as 
the surface brightness of most early-type galaxies. As a result, the 
surface brightness of the Abel models cannot accurately be fitted 
with smooth parameterizations such as the Multi-Gaussian Expan- 
sion method (MGE; Monnet et al. 1992; Emsellem et al. 1994), 
unless negative Gaussians and/or strong twists are invoked. This 
would not only be unnatural, but also incorrectly limit the possible 
deprojections. Nevertheless, the MGE parameterization allows for 
a straightforward deprojection and an efficient computation of the 
gravitational potential that correspond to the resulting (triaxial) in- 
trinsic density. Therefore, we enforce a less accurate but smooth 
MGE parameterizations of the surface brightness of the Abel models, 
and show in § 4.4 that this does not affect our main results. 
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3 CONSTRUCTING THE SCHWARZSCHILD MODELS 

Based on the simulated observations of each of the thirteen Abel 
models, we now construct triaxial Schwarzschild (1979) models 
using our numerical code described in detail in vdB08. The code re- 
constructs the intrinsic dynamical structure of a collisionless stellar 
system with an arbitrary triaxial geometry by finding the superposi- 
tion of orbits that best fits simultaneously the (observed) photometry 
and stellar kinematics in a self-consistent way. These models are 
significantly different from axisymmetric Schwarzschild models 
(e.g. van der Marel et al. 1998; Gebhardt et al. 2000; Valluri et al. 
2004; Thomas et al. 2004; Cappellari et al. 2006), as they allow for 
triaxial intrinsic shapes and the associated multiple distinct orbit 
families. Hence, they are applicable to many of the early-type galax- 
ies, including very boxy ones, and their triaxial dark haloes. The 
models presented in this paper use the default parameters as set in 
vdB08, including the number of orbits and the dimensions of the 
intrinsic polar mass grid. No regularisation was used throughout 
this paper, to ensure that the results unbiased towards any form of 
regularisation. 



3.1 Fitting the simulated observables of the Abel models 

The stellar kinematics fitted are the maps of the mean line-of-sight 
velocity V, velocity dispersion a and higher-order Gauss-Hermite 
velocity moments (up to ho) computed for each of the Abel models. 
In addition, we fit 2 the MGE parameterization of surface brightness 
distribution and the corresponding intrinsic luminosity density for 
a given viewing direction. Similarly, the intrinsic mass density and 
gravitational potential is inferred from the deprojected MGE pa- 
rameterization of the surface mass distribution. Because the surface 
mass distribution is not directly observable (except perhaps in case 
of gravitational lensing), we constructed the Abel models such the 
luminosity density mimics the mass density inside the observed 
region (see § 2 and vdV08). Essentially, we implicitly assume that 
mass follows light. Unlike the sometimes sharp transitions in the 
projected luminsity density as mentioned in § 2, the projected mass 
density of the Abel models is always smooth and hence accurately 
parameterized by a MGE. The accelerations in the (efficiently) com- 
puted MGE gravitational potential are accurate to < 1%, when 
compared to the input analytic Stackel potential. 

The types of orbits hosted by the gravitational potential depend 
strongly on the intrinsic shape of the galaxy. This means that the 
capability of the Schwarzschild models to fit the stellar kinemat- 
ics (and the photometry in a self-consistent way) also depends on 
the intrinsic shape (see for example Fig. 5 of vdB08). To obtain 
the intrinsic density distribution from deprojection of the MGE pa- 
rameterization three viewing angles are needed that are typically 
unknown: the polar ($) and azimuthal (ip) viewing angles, and the 
misalignment angle (ip) between the projected intrinsic short axis 
and the photometric minor axis. 

Searching through a uniform grid in these three angles is un- 
feasible since the deprojection is strongly non-linear. While in parts 
of the grid small changes in the viewing angles can result in strong 
variations in the intrinsic shape, in other regions the intrinsic shape 
remains nearly constant. Instead, we step through the possible de- 
projections by regularly sampling three intrinsic shape parameters: 
the axis ratios p and q, and a third parameter u that represents the 

2 Within the model the surface brightness and intrinsic density are actually 
treated as constraints instead of being directly fitted. 



scale-length compression factor. This has the advantage that we 
can sample the whole deprojection space with relatively few points 
and avoid sampling in parts of the viewing space where no changes 
occur, while increasing the sampling density where necessary (see 
also § 3.7 of vdB08). For the tests in this paper it is important to con- 
sider both the recovery of the viewing angles and the intrinsic shape 
separately. As a good constraint on the shape does not nessecarily 
mean a good constraint on the viewing angles and vice versa. 

The (p, q, u) parameters are not equivalent to the Stackel ps 
and qs from § 2: both describe the intrinsic shape, but the latter 
describes the intrinsic shape measured directly from the moments 
of inertia in the isochrone Stackel potential (see § 4. 1 in vdV08), 
whereas (p, q, u) is defined within the MGE system. 

The range of deprojections, and thus the sampling in the in- 
trinsic shape parameters, is different for each of the galaxy models 
presented in this paper, as they have a different surface brightness. 
Typically we have 125 deprojections in steps of 0.05 in (p, q, it), 
and for each deprojection we sample the mass-to-light ratio M/L 
at eleven different values. In this way, we construct around 1400 
Schwarzschild models to the simulated observables of each Abel 
model. Computing the orbit library - which has to be done for every 
deprojection - takes approximately six hours on a current single- 
core CPU. For each deprojection and M/L the fitting takes another 
hour. This means in total about 2000 CPU hours to search the whole 
parameter space per Abel model. 

3.2 Confidence criterion 

To determine the best-fit model parameters and corresponding uncer- 
tainties we compare for each Schwarzschild model its x 2 difference 
between the simulated observables and the corresponding model 
predictions, weighted with the errors in the observables. The (global) 
minimum in the x' 2 yields the best estimate of the model parameters. 
It is common practice, e.g., when dynamical models are used to 
measure the mass of a black hole, to assign error bars based on 
Ax 2 = X 2 ~ m i n (x 2 ) an d the number of free parameters M of 
the model. In our case we have M — 4 parameters (three view- 
ing angles and the mass-to-light ratio) after marginalizing over the 
orbital weights. And thus our assiociated formal errors become 
Ax 2 = 4.72 for 68% confidence (la). 

Whereas for the black hole mass measurements, with typically 
only a few observations inside the sphere of influence of the black 
hole, this works well, a different criterion is needed in our case of 
measuring the intrinsic shape based on (simulated) integral-field 
spectroscopy observations. The integral-field spectroscopy obser- 
vations that are used typically consist of several hundred (or more) 
spectra. From each spectrum four (or more) velocity moments are ex- 
tracted, yielding a large number of individual constraints N > 10 4 , 
which are fitted by the Schwarzschild model. Under these circum- 
stances the uncertainties in the x 2 itself become important. The 
X 2 itself has a standard deviation of \f2(N — M), which becomes 
non-neglibible in our case. Moreover, we expect that systematic 
errors in the observations (e.g. due to stellar template mismatch) can 
cause deviations in the x of the order of (or more than) the formal 
la. 

Instead, we derive uncertainties on the best-fit parameters based 
on the expected standard deviation of x 2 , which we approximate to 
be y/2N because N M. Preliminary validation of this approach 
is shown in vdB08, where we construct triaxial Schwarzschild mod- 
els of one Abel model (similar to ST4) and the elliptical galaxy 
NGC 4365 with an (apparent) kinematically decoupled core. In this 
case and in the current study, the viewing angles of the Abel models 
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are known, so that this criterion for the standard error bars can be 
validated. 



4 RESULTS OF THE DYNAMICAL MODELS 

In Figures 1-5, we present for each Abel model the simulated ob- 
servables, the best-fit Schwarzschild model, and the Schwarzschild 
model closest to the input parameters. The difference between the 
latter two Schwarzschild models reveals the (often small) differ- 
ences due to the typical uncertainty allowed in the fitted stellar 
kinematics (e.g. template mismatch errors). The contour plots on the 
right indicate the confidence on the recovered M/L, the shape and 
the corresponding viewing angles. Besides the best-fit model (red 
crosses), they also show the location of the true input values (green 
diamonds). The white areas do not have a feasible deprojection. We 
first discuss the individual cases. 

4.1 The oblate models 

For all six oblate Abel models (SOI through F03) shown in Fig- 
ures 1 and 2, the Schwarzschild models around the minimum in \ 2 
are (very) good fits with reduced x 2 P er degrees-of-freedom close to 
unity. The shape is recovered well within the la confidence region. 
The edge-on and intermediate models (SOI, S02, FOl, F02) have 
a strong constraint on p > 0.95, fully excluding strongly triaxial 
geometries, while q is generally unconstrained, apart from the lower 
limit q > 1 — e implied by the largest value of the observed ellip- 
ticity e. Hence, we cannot constrain the inclination of these (near) 
axisymmetric systems. The differences in the fitted kinematics be- 
tween the best-fitting model and the model with the correct input 
parameters (see Fig. 1) are so small that they are negliable compared 
to possible (systematic) errors in real observations. 

As expected, the face-on slow-rotator model S03 has the worst 
recovery with no significant constraint on its shape and a wide range 
in M/L allowed. All the other oblate models have a constraint on 
M/L that is consistent with the input value. The constraints on all 
the fast-rotator oblate models are tighter than those on the slow- 
rotator oblate models because the mean velocity amplitude in the 
fast-rotator models is larger, which effectively shrinks the simulated 
observational errors. 

For all oblate models the flattening q of the best-fit model is 
never far removed from the input value go: Aq = \q — qo\ < 0.08. 
Even though the viewing direction, or the inclination in this case, 
is mostly unconstrained, the latter indicates that the intrinsic shape 
is not fully degenerate. Observations that are more accurate and/or 
span a larger field-of-view are expected to yield tighter constraints on 
the intrinsic shape. When we fitted the surface brightness of the Abel 
models directly instead of the overly smooth MGE parameterization, 
the constraints on the shape improve, but only mildly, showing that 
our MGE parameterization is sufficient for the oblate models. 

4.2 The prolate models 

The results for the (slow-rotator) prolate models (SP1, SP2, SP3) 
are shown in Fig. 3. The constraint on the shape is better than for 
the oblate models and the true shapes and M/L lie within or on the 
la confidence contour. In all cases the axis ratios are well recovered 
with allowed ranges Ap < 0.07 and Aq < 0.07. Even the end-on 
model SP3 has a good constraint on the shape, excluding more than 
75% of the allowed range in p and q. The constraint on the shape for 
prolate systems is better than for the oblate systems due to the fact 



that there are less (near) prolate deprojections possible than oblate 
ones. Similarly, the limits on the inclination from the photometry are 
more restrictive in the prolate than oblate case, but still the regular 
kinematics do not really further constrain the inclination. 

4.3 The triaxial models 

Model ST1 (Fig. 1) is a non-rotating round galaxy seen at an inter- 
mediate viewing angle. Like the face-on slow-rotator oblate model 
S03 - which shows almost no rotation - the triaxial Schwarzschild 
models are capable of reproducing all features in the simulated ob- 
servables at nearly all viewing angles. The M/L is recovered with a 
similar accuracy as in the case of the oblate models, despite the lack 
of constraint on the intrinsic shape. Even so, the best-fit p and q are 
with 10 per cent from the input values, while the formally allowed 
ranges are larger with Ap ~ 0.30 and Aq ~ 0.30. 

The zero velocity curve of the stellar kinematics in model ST2 
seen in Fig. 4 is twisted with respect to the photometric minor 
axis. Since this is only possible in a non-axisymmetric system, we 
might expect a rather good recovery of the intrinsic shape. However, 
the constraints are not very good with ranges Ap ~ 0.15 and 
Aq ~ 0.25, and comparable with those on the shape of the slow 
rotator oblate models SOI and S02. Model ST3 (Fig. 4) has a 
central rotating component and no rotation in the outer parts. The 
resulting constraint on the intrinsic shape with ranges Ap ~ 0.15 
and Aq ~ 0.20 is only slightly better than model ST2. However, as 
we show below in § 4.4 and in Fig. 6, taking into account the strong 
photometric twist improves the recovery of the intrinsic shape as 
well as the (weak) constraint on the viewing angles. 

Model ST4 (Fig. 5) contains a central KDC like model ST3, and 
additionally shows significant rotation in the outer parts. Although 
the inner part of the surface brightness is not well reproduced by 
the smooth MGE parameterization, unlike model ST3, the intrinsic 
shape is nevertheless accurately recovered and constrained, with 
ranges Ap < 0.05 and Aq < 0.05. The triaxial Abel model used 
in vdB08 and vdV08 is similar to model ST4, but at a different 
viewing direction and with more photometric twist in the MGE 
parameterization, resulting in (slightly) tighter constraints on the 
intrinsic shape as well as on the viewing angles. 

4.4 Fitting the non-smooth surface brightness directly 

As mentioned in § 2, the surface brightness of the Abel models can 
show a sudden transition between DF components, which cannot be 
accurately described by a smooth MGE parameterization without 
invoking negative Gaussians and/or strong photometric twists. A 
way to avoid this is to bypass the MGE and use the observed sur- 
face brightness directly. This means that the intrinsic density is not 
constrained anymore, since it follows from deprojecting the MGE 
parameterization. In Fig. 6, we present the results for model ST3, 
which shows most evidently both a central depression as well as a 
strong photometric twist in its surface brightness. The intrinsic shape 
is still correctly recovered, but now the best-fit model is much closer 
to the input values. Also many of the deprojections with p > 0.93 or 
q > 0.67 are excluded, leading to a significantly tighter constraint 
on the intrinsic shape then when the MGE parameterization is used. 
In Fig. 6 we also show the same for ST2, for which the effect is 
less pronounced but still significant: the best-fit model is now much 
closer with p and q less than 0.05 away from their input values, but 
the confidence intervals are unchanged. 

This shows that the MGE representation we used in § 4.3 is too 
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Figure 1. Dynamical modelling of fast rotator oblate models FOI, F02 and F03. Each plot shows the stellar kinematics of the Abel model (top), the best-fit 
Schwarschild model (middle), and the fit at the input parameters (bottom). The contours on the surface brightness image (leftmost panel) show the observed 
surface brightness (top), the surface brightness of the MGE model (bottom) and the difference between the two (middle). Colums 2 through 5 show the 
line-of-sight velocity moments: V (-110. . . 110 km s -1 ), a (100. . . 260 km s -1 ), /13 and /14 (-0.1. . . 0. 1). The panels on the right side show marginalized 
confidence regions of the viewing angles j9 and ip (in Lambert projection seen down the 2-axis, so that # runs radially and (p runs azimuthally), p versus q, and p 
versus M/L. The inner red and outer black contour signify lc and 3cr confidence, repectively. The red cross and green diamond highlight the best-fit parameters 
and the original input parameters, respectively. 
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Figure 4. Slow rotator triaxial models STI, ST2 and ST3. See Fig. 1 for description. 
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Figure 6. The shape recovery of model ST2 (left) and ST3 (right), when 
the observed surface brightness is fitted directly, bypassing the (too) smooth 
MGE parameterization. The presentation of the models here is identical to 
the ones in in Fig. 4. 

smooth, adding irrelevant deprojections that (artificially) increase 
the error bars on p and q. However, there is a delicate balance 
between reproducing too little and too much photometric twist. 
Fitting unrealistic sharp features (due to for example dust in real 
observations), may exclude allowed deprojections or even the true 
underlying intrinsic shape. To avoid the latter, the fit to the surface 
brightness should be as round as possible, while still accurately 
describing the observed photometry. 



5 INTRINSIC VELOCITY MOMENTS 

In vdV08, we have shown that given the correct input viewing direc- 
tion, the internal dynamics and three-integral distribution function 
can be recovered. Here, we investigate what happens to the intrinsic 



velocity moments when the viewing angles are allowed to deviate 
from the input values. We expect the intrinsic velocity moments, 
in terms of V/a and velocity anisotropy, to vary when the view- 
ing direction and hence intrinsic shape changes. Still, if we stay 
within the la confidence regions of the fitted intrinsic shape pa- 
rameters, these variations might be limited due to the constraints 
coming from the two-dimensional kinematics. We use model ST3 
from the previous section as an illustration. In Fig. 7, we show the 
intrinsic velocity moments at the input values (p, q) = (0.85, 0.66) 
and at (p, q) = (0.90, 0.50), which is at the (lower) boundary of 
the corresponding la confidence region (see also Fig. 6). While the 
mass density of the two models (indicated by the black contours) 
is clearly different, the intrinsic velocity moments are quite similar. 
This suggests that as long as the intrinsic shape is reasonably well 
constrained, the internal dynamics are also reliably recovered. 



6 DISCUSSION AND CONCLUSIONS 

We constructed thirteen galaxy models with simulated photom- 
etry and two-dimensional stellar kinematics that represent the 
SAURON observations of early-type galaxies, from oblate fast rota- 
tor to triaxial slow rotator. We fitted realistic three-integral triaxial 
Schwarzschild models to the simulated observables of these Abel 
models and we used the x 2 difference to investigate how well the 
intrinsic shape as well as the mass-to-light ratio can be recovered. 
Due to the large number of observational constraints (typically 
N > 10 4 ) with possible systematic errors (e.g. stellar template 
mismatch) in the case of real integral-field observations, the com- 
monly used A\ 2 — X 2 ~ mm (x 2 ) ls not useful to estimate the 
uncertainties on the model parameters. Instead, we showed that the 
standard deviation %/2N works well as a confidence criterion. 

The simulated observations of the (near) oblate Abel models 
that are not viewed close to face-on strongly excluded triaxial in- 
trinsic shapes with 1 — p < 0.05. The flattening q is almost only 
constrained by the photometry, although the models with more pro- 
jected rotation provide a slightly tighter constraint since the relative 
errors in the odd velocity moments are smaller. This means that for 
fast-rotator early-type galaxies, which seem to be consistent with 
oblate axisymmetry, the inclination is not significantly further con- 
strained by current two-dimensional stellar kinematics, consistent 
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Figure 7. The intrinsic velocity moments of triaxial Schwarzschild models fitted to the simulated observables of the Abel model ST3. In the left panel, the 
intrinsic shape parameters are at the true input values, while in the right panel they are at the boundary of the corresponding lcr confidence region. The colour 
represents the (intrinsic) V/cr, and the ellipses are cross-sections of the velocity ellipsoid with the (x, z) plane, while the crosses indicate the (relative) size of 
the velocity ellipsoid in the perpendicular (j/-axis) direction (see § 5.3 of vdV08 for further details). The black curves are contours of constant mass density in 
steps of 1 mag. 



with Krajnovic et aL (2005). Similarly, for the (near) prolate Abel 
models that are not viewed close to end-on a triaxial intrinsic shape 
is strongly excluded with p — q < 0.05, but the viewing direction is 
only weakly constrained. 

Unlike axisymmetric galaxies, a triaxial galaxy can show twists 
in both the surface brightness and velocity field (and in higher-order 
velocity moments), as well as misalignment between the kinematic 
(rotation) axis and the photometric minor axis. The stronger such 
features are in the observations, the more accurate the intrinsic 
triaxial shape can be recovered, often better than in the axisymmetric 
case. In particular, for the triaxial Abel models with a central KDC, 
p and q can be recovered within 5%, especially when (misaligned) 
rotation is present in the outer parts. Since such kpc-size KDCs 
seem typical for the slow-rotator early-type galaxies, we expect that 
their intrinsic shape can be well estimated from dynamical modeling. 
However, when inferring the mass model from a fit to the surface 
brightness, the fit should be as round as possible to not exclude 
allowed deprojections, but still accurately describe the photometry, 
including twists, to avoid artificially increasing the range of allowed 
intrinsic shapes. 

Near face-on oblate galaxies and end-on prolate galaxies ap- 
pear round on the sky and show almost no (projected) rotation. We 
showed that in these cases, as well as for triaxial galaxies with no 
intrinsic (net) rotation, the best-fit p and q are still remarkably close 
to their true values (typically within 10%), but the formal uncer- 
tainties are much larger (typically around 30%). This means that 
for galaxies that appear round on the sky and show little rotation, 
the only hope of recovering the intrinsic shape and corresponding 
viewing direction is through additional constraints such as provided 
by a thin embedded disc (in either dust, gas or stars). 

In all cases, the M/L is correctly recovered, within about 10% 
when the intrinsic shape is well recovered, increasing up to 20% 
for the models with a poor constraint on the intrinsic shape. This 
is in line with Thomas et al. (2007), who showed that applying 



axisymmetric Schwarzschild models to triaxial merger remnants 
results in a larger scatter in the recovered M/L compared to oblate 
remnants. However, as long as generic dynamical models such as 
our triaxial Schwarzschild models are used, the recovery of the 
M/L is nearly independent of the underlying intrinsic shape. This is 
not unexpected given the good correlation between dynamical and 
virial M/L estimates (e.g. Cappellari et al. 2006), and the observed 
tight scaling relations in early-type galaxies such as the fundamental 
plane (e.g. Dressier et al. 1987; Djorgovski & Davis 1987). 

In real galaxies the photometry and the kinematics are usually 
measured independently and as a result their relative alignment is 
not perfectly known. This can create a small kinematic misalignment 
in an otherwise perfectly axisymmetric galaxy, or vice-versa. To 
test the significance of this we re-ran the fast rotator oblate model 
F02 with a positive and negative misalignment of 2°, emulating a 
typical measurement error on the PA (in a fairly round galaxy). The 
resulting recovered parameters and their errors were very similar 
to the original F02 model, with the \ 2 om y marginally higher 
and the 3u confidence region slightly extended. This shows that 
the uncertainty in the relative PA in these observations does not 
effect the recovery of the intrinsic shape nor the recovery of the 
mass-to-light ratio. 

Furthermore, we ran seperate models of F02 and ST4 in which 
we doubled the number of orbits. The resulting models, including 
the \ 2 contours and the recovered kinematics, where nearly identical 
showing that the number of orbits used in our models is sufficient. 

We conclude that for (near) axisymmetric galaxies the com- 
bination of photometric and two-dimensional stellar kinematic ob- 
servations can strongly exclude triaxiality, but regular kinematics 
do not further tighten the intrinsic flattening significantly. As a re- 
sult, we expect the inclination of fast-rotator early-type galaxies to 
remain nearly unconstrained above the photometric lower limit, al- 
though better observations do seem to slightly decrease the range in 
inclinations (§4.1). Triaxial galaxies can have additional complexity 
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in both the observed photometry and kinematics, such as the central 
KDCs observed in many slow-rotator early-type galaxies, which 
allows the intrinsic shape to be accurately recovered. The intrinsic 
shape of round galaxies with no significant rotation is degenerate, 
unless additional constraints such as from a thin disk are available. 
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